function plotTermPremia(outEstim)

mat = 40; %Maturity for Term Premia
setupFilter     = outEstim.setupFilter;
%timeIndexPlot   = 1960:10:2020;  
dateTimes       = outEstim.setupFilter.timeIndex;
timeIndexPlot   = round(year(dateTimes(1,:))/10)*10:5:round(year(dateTimes(end,:))/10)*10;

fontSizeValue   = 14;
figureSize      =  [0 0 1300 600+300]; % (0,0) size in x size in y

termPrem        = outEstim.termPrem;
termPremReal    = outEstim.termPremReal;

%% Plot Nominal Term premium 
version = 1;
figure('Name','Nominal Term Premium','NumberTitle','off','Position',figureSize)
hold on
NBER        = outEstim.setupFilter.NBERrec(setupFilter.numInitial+1:end);
yMax        = max([termPrem.TP(mat,setupFilter.numInitial+1:end),setupFilter.acm_tp.tp10y']);
yMin        = min([termPrem.TP(mat,setupFilter.numInitial+1:end),setupFilter.acm_tp.tp10y']);
dates       = setupFilter.timeIndex(setupFilter.numInitial+1:end);
if version == 1
    xaxis       = setupFilter.timeIndex(setupFilter.numInitial+1:end);
    nn          = length(dates);
    idx         = nn+setupFilter.numInitial:-1:setupFilter.numInitial+1;
    reverseaxis = setupFilter.timeIndex(idx);
    index       = 1:nn;
    reversIndex = nn:-1:1;
    grpyat = [xaxis       100*yMax*NBER(index,1);...
        reverseaxis NBER(reversIndex,1)*yMin];
    patch(grpyat(:,1),grpyat(:,2),[0.8 0.8 0.8],'edgecolor',[0.6 0.6 0.6]);    
else
    %For shading above zero
    X  = [dates',fliplr(dates)'];
    Y1 = [100*yMax*NBER' 0*NBER'];
    p2 = fill(X,Y1,[0.8 0.8 0.8]);
    set(p2,'edgecolor',[0.8 0.8 0.8]);
    
    % For shading below zero
    Y2 = [100*yMin*NBER' 0*NBER'];
    p3 = fill(X,Y2,[0.8 0.8 0.8]);
    set(p3,'edgecolor',[0.8 0.8 0.8]);
end
hold on
tmp1 = plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),termPrem.TP(mat,setupFilter.numInitial+1:end),'-k','LineWidth',2);
tmp2 = plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),...
    setupFilter.acm_tp.tp10y(setupFilter.numInitial+1:end),'-k');
legend([tmp1(1), tmp2(1)],{'$M^{M,CS}$','Adrian et al (2013)'},'Orientation','horizontal','Location','NorthEast','EdgeColor',[1 1 1],...
               'interpreter','Latex','FontSize',fontSizeValue+4);
hold off

%axis tight
axis([dates(1) dates(end) yMin yMax])

datetick('x','yyyy','keeplimits','keepticks')
xticks(datenum(datetime(timeIndexPlot,1,1)))
xticklabels(num2cell(timeIndexPlot))
set(gca,'FontSize',fontSizeValue+2)
corrTPnom_ACM = corr(termPrem.TP(mat,setupFilter.numInitial+1:end)',setupFilter.acm_tp.tp10y(setupFilter.numInitial+1:end))
%title([num2str(mat/4),'-year term premium: Mean = ',num2str(termPrem.MeanTP(mat,1),3)],'FontSize',fontSizeValue+2)
%title([num2str(mat/4),'-year term premium '],'FontSize',fontSizeValue+2)

%% Plot Real Term premium 
version = 1;
figure('Name','Real Term Premium','NumberTitle','off','Position',figureSize)
hold on
yMax        = max(termPremReal.TP(mat,setupFilter.numInitial+1:end));
yMin        = min(termPremReal.TP(mat,setupFilter.numInitial+1:end));
if version == 1
    xaxis       = setupFilter.timeIndex(setupFilter.numInitial+1:end);
    nn          = length(dates);
    idx         = nn+setupFilter.numInitial:-1:setupFilter.numInitial+1;
    reverseaxis = setupFilter.timeIndex(idx);
    index       = 1:nn;
    reversIndex = nn:-1:1;
    grpyat = [xaxis       100*yMax*NBER(index,1);...
        reverseaxis NBER(reversIndex,1)*yMin];
    patch(grpyat(:,1),grpyat(:,2),[0.8 0.8 0.8],'edgecolor',[0.6 0.6 0.6]);    
else
    %For shading above zero
    X  = [dates',fliplr(dates)'];
    Y1 = [100*yMax*NBER' 0*NBER'];
    p2 = fill(X,Y1,[0.8 0.8 0.8]);
    set(p2,'edgecolor',[0.8 0.8 0.8]);
    
    % For shading below zero
    Y2 = [100*yMin*NBER' 0*NBER'];
    p3 = fill(X,Y2,[0.8 0.8 0.8]);
    set(p3,'edgecolor',[0.8 0.8 0.8]);
end
plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),termPremReal.TP(mat,setupFilter.numInitial+1:end),'-k')

%axis tight
axis([dates(1) dates(end) yMin yMax])
datetick('x','yyyy','keeplimits','keepticks')
xticks(datenum(datetime(timeIndexPlot,1,1)))
xticklabels(num2cell(timeIndexPlot))
set(gca,'FontSize',fontSizeValue)
title([num2str(mat/4),'-year real term premium: Mean = ',num2str(termPremReal.MeanTP(mat,1),3)],...
    'FontSize',fontSizeValue+2)

%% Plot Inflation Risk Term premium 
version = 1;
figure('Name','Inflation Risk Term Premium','NumberTitle','off','Position',figureSize)
hold on
inflTP      = outEstim.termPrem.TP - outEstim.termPremReal.TP;
yMax        = max(inflTP(mat,setupFilter.numInitial+1:end));
yMin        = min(inflTP(mat,setupFilter.numInitial+1:end));
if version == 1
    xaxis       = setupFilter.timeIndex(setupFilter.numInitial+1:end);
    nn          = length(dates);
    idx         = nn+setupFilter.numInitial:-1:setupFilter.numInitial+1;
    reverseaxis = setupFilter.timeIndex(idx);
    index       = 1:nn;
    reversIndex = nn:-1:1;
    grpyat = [xaxis       100*yMax*NBER(index,1);...
        reverseaxis NBER(reversIndex,1)*yMin];
    patch(grpyat(:,1),grpyat(:,2),[0.8 0.8 0.8],'edgecolor',[0.6 0.6 0.6]);    
else
    %For shading above zero
    X  = [dates',fliplr(dates)'];
    Y1 = [100*yMax*NBER' 0*NBER'];
    p2 = fill(X,Y1,[0.8 0.8 0.8]);
    set(p2,'edgecolor',[0.8 0.8 0.8]);

    % For shading below zero
    Y2 = [100*yMin*NBER' 0*NBER'];
    p3 = fill(X,Y2,[0.8 0.8 0.8]);
    set(p3,'edgecolor',[0.8 0.8 0.8]);
end
plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),inflTP(mat,setupFilter.numInitial+1:end),'-k')

%axis tight
axis([dates(1) dates(end) yMin yMax])
datetick('x','yyyy','keeplimits','keepticks')
xticks(datenum(datetime(timeIndexPlot,1,1)))
xticklabels(num2cell(timeIndexPlot))
set(gca,'FontSize',fontSizeValue)
title([num2str(mat/4),'-year inflation risk term premium: Mean = ',num2str(mean(inflTP(mat,:)),3),...
' Std = ',num2str(std(inflTP(mat,:)),3)],'FontSize',fontSizeValue+2)


%% Shock Decomposition of Nominal Term Premia
timeIndexPlot   = round(year(dateTimes(1,:))/10)*10:10:round(year(dateTimes(end,:))/10)*10;
version = 1;
figure('Name','Nominal Term Premium: Decomposition','NumberTitle','off','Position',figureSize)
hold on
NBER        = outEstim.setupFilter.NBERrec(setupFilter.numInitial+1:end);
yMax        = max(max(termPrem.TPdecom(mat,setupFilter.numInitial+1:end,:)));
yMin        = min(min(termPrem.TPdecom(mat,setupFilter.numInitial+1:end,:)));
dates       = setupFilter.timeIndex(setupFilter.numInitial+1:end);

termPremDeMean_TPdecom = termPrem.TPdecom - repmat(mean(termPrem.TPdecom,2),1,size(termPrem.TPdecom,2),1);
totalAbs      = sum(abs(termPremDeMean_TPdecom(mat,setupFilter.numInitial+1:end,:)),3);
fracShockAbs  = squeeze(sum(abs(termPremDeMean_TPdecom(mat,setupFilter.numInitial+1:end,:)))./sum(totalAbs));
totalRMSE     = sum((termPremDeMean_TPdecom(mat,setupFilter.numInitial+1:end,:).^2),3);
fracShockMSE  = squeeze(sum(termPremDeMean_TPdecom(mat,setupFilter.numInitial+1:end,:).^2)./sum(totalRMSE));

for i=1:outEstim.model.ne
    subplot(2,3,i)
    if version == 1
        xaxis       = setupFilter.timeIndex(setupFilter.numInitial+1:end);
        nn          = length(dates);
        idx         = nn+setupFilter.numInitial:-1:setupFilter.numInitial+1;
        reverseaxis = setupFilter.timeIndex(idx);
        index       = 1:nn;
        reversIndex = nn:-1:1;
        grpyat = [xaxis       100*yMax*NBER(index,1);...
            reverseaxis NBER(reversIndex,1)*yMin*100];
        patch(grpyat(:,1),grpyat(:,2),[0.8 0.8 0.8],'edgecolor',[0.6 0.6 0.6]);
    else
        %For shading above zero
        X  = [dates',fliplr(dates)'];
        Y1 = [100*yMax*NBER' 0*NBER'];
        p2 = fill(X,Y1,[0.8 0.8 0.8]);
        set(p2,'edgecolor',[0.8 0.8 0.8]);
        
        % For shading below zero
        Y2 = [100*yMin*NBER' 0*NBER'];
        p3 = fill(X,Y2,[0.8 0.8 0.8]);
        set(p3,'edgecolor',[0.8 0.8 0.8]);
    end

    hold on
    %plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),...
    %    termPrem.TP(mat,setupFilter.numInitial+1:end)-termPrem.TPdecom(mat,setupFilter.numInitial+1:end,i),'-k');
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),...
        termPrem.TPdecom(mat,setupFilter.numInitial+1:end,i),'-k');
    hold off
    axis([dates(1) dates(end) yMin yMax])
    datetick('x','yyyy','keeplimits','keepticks')
    xticks(datenum(datetime(timeIndexPlot,1,1)))
    xticklabels(num2cell(timeIndexPlot))
    set(gca,'FontSize',fontSizeValue)
    title([outEstim.model.labelx{outEstim.model.nx1+i},...
        ': pct($\Delta TP^i_t$) = ',num2str(round(100*fracShockAbs(i),1))],...
        'Interpreter','Latex','FontSize',fontSizeValue+2)
    %title([outEstim.model.labelne{i},': pct($\Delta TP_t$) = ',num2str(round(100*fracShockRMSE(i),1))],'Interpreter','Latex','FontSize',fontSizeValue)
end


end